Finding the structure of phosphorus in the phase IV 
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We have explored the unknown structure of the phosphorus in the phase IV (P-IV phase) based 
on the first-principles calculations using the metadynamics simulation method. Starting from the 
simple cubic structure, we found a new modulated structure of a monoclinic lattice. The modulation 
is crucial to the stability of the structure. Refining further the structure by changing the modulation 
period, we have found the structure which shows the X-ray powder pattern in the best agreement 
with the experimental one. We can not exclude the possibility that the unknown structure of the 
phase IV of phosphorus is an incommensurately modulated one. 
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Recent progress in high-pressure physics strengthened 
our recognition on the variety of structures of materi- 
als. Unexpectedly interesting structures were found in 
the high-pressure experiments. The modulated struc- 
ture is one of curious structures often found in the high- 
pressure phases of elements. Improvement of the high- 
pressure techniques has also made it possible to specify 
the structure stabilized only in a narrovif pressure range. 
The modulation of a crystal was found in the group Vb 
elements including As, Sb and Bi 0, and in group 
VIb elements including S 0, Se ||] and Te 0. Modu- 
lated structures were also reported in halogens, I and 
Br 0. 

Non availability of sufficient experimental data under 
very high pressure constrains high pressure study. There, 
people often encounter difficulties in determination of the 
lattice structure only from the experimental data. The 
theoretical approach gives alternative information on the 
same problem, where accuracy in the determination of a 
crystal structure is known to be enough, if one utilizes the 
first-principles calculation. However, limitation in com- 
putational resources sometimes prohibits us to perform 
full search of the true structure. 

The phase IV of the phosphorus (P-IV) is one of the ex- 
amples in which the structure has not been determined. 
Observation of the phase was first reported by Akahama 
et al. in 1999. They showed an appearance of the sim- 
ple hexagonal (sh) phase, i.e. the phase V (P-V), which 
is stabilized above 137GPa. The third phase (P-III) ap- 
pearing in the sequence of pressure-induced transforma- 
tions at low temperatures is the simple cubic (sc) phase. 
As an intermediate phase between sc and sh, the P-IV 
was detected in the X-ray diffraction data. Experimen- 
tally, however, the structure has not been identified. The 
ordinary Rietveld analysis starting from the knowledge 
only of the monoclinic symmetry has not been successful 
owing to a probable complexity of the lattice. The bcc 
structure (phase VI) is theoretically predicted and 
identified at even higher pressure 262GPa We need 
a guess of the lattice structure or a pseudo crystal. 



Several structures have been tested as candidates for 
the P-IV. Ahuja considered a structure of space group 
Imma 11]. Ehlers and Christensen studied relative sta- 
bility of Ba-IV structure of P, which is a kind of mod- 
ulated structure, in the pressure range from 100 to 200 
GPa|l^. In spite of these intensive studies, the structure 
of the P-IV has not been identified. 

To explore the P-IV, we adopted the following strat- 
egy in our theoretical study and used the metadynamics 
simulation in the first-principles calculations. This trial 
was done with a relatively small simulation cell to reduce 
the computational time. The simulation was planned, 
however, to be able to detect possible signal of the struc- 
tural phase transformation. For the obtained structure, 
we checked the relative stability against the sc and sh 
phases. Next we considered some model structures to 
find more refined structure. The structural optimization 
was done for each model structures. The calculated X-ray 
powder patterns of the optimum structure are compared 
with that of the experimental one. 

In the study of the metadynamics simulation, which 
was first introduced by Laio and Parrinello 0, Il3 |. 
we use the Gibbs free energy (GFE) depending upon 
a shape of a simulation cell. Following the prescrip- 
tion by Martohak HI, we consider the GFE, G(h*) = 
6*0(11*) -I- Gg{h*), where G'g(h') is the artificial potential 
defined by the following Gaussian-type function. 
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The superscripts t and t' denote the current meta-step 
and the previous one, respectively. The quantities W and 
Sh represent the weight and the width of the Gaussian- 
type function, respectively. The matrix h is defined by 
the vectors defining the simulation cell and h = (a, 6, c) , 
where a, b and c are lattice vectors. To eliminate the free 
rotation of the system, only the symmetric part of the 
matrix h is updated, which reduces the number of the 
independent variables to 6. 

The update of the matrix h is made by the steepest de- 
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scent method using the driving force F and regarding 5h 
as stepping parameter. The driving force F is obtained 
as the sum of the original driving force F^ = —dGo/dh, 
and the Gaussian driving force Fg — —dGg/dh. The 
force Fo can be expressed by an internal pressure tensor 
p, external pressure P, and the matrix h 0|. One step 
of updating h is defined as one meta-step. 

At each meta-step, in order to equilibrate the system 
and to estimate p, the conventional molecular dynam- 
ics (MD) simulations have to be done with the shape of 
the simulation cell fixed. The internal pressure tensor 
and the atomic positions at each meta-step can be taken 
from the output of any constant-pressure MD codes of 
the first-principles calculation [isj . 

The above artificial potential means that if the current 
h has been visited time after time, which occurs when the 
system is fluctuating around the local minimum of the 
GFE surface, Gg accumulates to a large value and the 
well is gradually filled with the artificial potential Gg . 

For the simulation of phosphorus, we used the den- 
sity functional theory in a local density approximation 
and a norm-conserving pscudopotential, where we em- 
ployed the expression of Perdew and Zunger for the 
exchange and correlation energy functional. We checked 
the pscudopotential comparing the calculated equation 
of states with the experiments in other phases |3, M, USM ■ 
We started with the cubic simulation cell whose edge was 
4.26 A and 8 phosphorus atoms are set at the positions 
which make the sc lattice. We performed the k-space 
integration using 8x8x8 mesh points in the first Bril- 
louin zone and set the energy cut-off of the plane wave 
basis at 40 Ry. We set the external pressure at 120 GPa 
in the conventional constant-pressure MD, since the P- 
IV is observed at this pressure. In order to equilibrate 
the system, we ran this MD simulation for 200 steps at 
each meta-step and we calculated the average internal 
pressure tensor from the latter half of 100 steps. In order 
to perform this metadynamics simulation, we used the 
cluster machines. 

FigureQlshows the evolution of the cell volume and the 
three angles among the lattice vectors of the simulation 
cell. First we ran the simulations with the Gaussian po- 
tential off, by setting W = mRy, for 38 meta-steps in 
order to check whether the sc structure is stable or not. 
During these initial meta-steps, the values of the angles 
and the volume were nearly maintained at those of the 
sc lattice. This means that the sc structure resides on a 
local minimum and it is separated from the other local 
minima by some barriers. We note that the metadynam- 
ics simulation with the Gaussian-type potential switched 
off is nearly equivalent to the conventional variable-cell 
MD simulations. 

In order to explore the metastable structures beyond 
the potential barrier we switched the Gaussian-type po- 
tential on at 39th meta-step with of 1 mRy and Sh 
of 20 mA. As a result, one of the three angles started 
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FIG. 1: Evolution of the simulation cell volume and the 
angles among lattice vectors. The Gaussian-type potential 
was switched on at 39th meta-step and off again after 71st 
meta-step. 
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FIG. 2: A structure obtained by the first-principles metady- 
namics simulation, a, b and care lattice vectors of the simula- 
tion cell, and a = 4.22 A, 6 = 4.15 A, c = 4.22 A, a = 90.86°, 
P = 97.76° and 7 = 90.26°. When we look at the ac plane 
along the 6-axis, the planes are displaced alternatively in the 
direction [110]. 



to increase after around 50th meta-step and the volume 
began to decrease dramatically. After those changes, we 
switched the Gaussian-type potential off at 71st meta- 
steps again in order to check whether the system had al- 
ready surmounted the barrier and had moved to a neigh- 
boring local minimum. If the system had not crossed 
the barrier yet, the angle and the volume would have re- 
turned to the starting values of the sc lattice, which are 
approximately 90° and 74 A'^, respectively. After 71th 
meta-step, however, the volume fluctuated around 73 A^ 
and three angles also continued to fluctuate around 90°, 
98° and 90°. This behavior shows that the sc structure 
transformed into another metastable one. 

Figure [5] shows the structure obtained by the above 
run, where a, b and c are the lattice vectors of the sim- 
ulation cell. This structure has a unit cell with the lat- 
tice parameters a = 4.22 A, = 4.15 A, c = 4.22 A, 
a = 90.86°, P = 97.76° and 7 = 90.26°. The left hand 
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FIG. 3: Modulated structures with diflterent modulation pe- 
riods along the fe-axis. The pattern Ml is of a non-modulated 
structure, and M2 the structure obtained by our metadynam- 
ics simulation (ABAB---). The modulation periods of M4 
(ABAC- ■ ■) and M8 (ABCBADED- ■ ■) patterns are twice and 
four times as long as that of the M2, respectively. 



TABLE I: Comparison of the total energies per atom for 
a non-modulated structure (Ml) and three modulated struc- 
tures (M2, M4 and M8). We used the unit ceU with the lattice 
vectors, a/2, 46 and c/2 in the calculation of the total energy 
for all of these four structures. 



Structure 


Total energy [Ry/atom] 


Ml 


-13.1624 


M2 


-13.1638 


M4 


-13.1642 


M8 


-13.1637 
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FIG. 4: Comparison of the experimentally obtained X-ray 
powder pattern of the P-IV Q with those of the Ml, M2, M4 
and MS modulated structures. Theoretical powder patterns of 
our structures are obtained by the use of RIETAN-2000 [T^ . 
In the experimental pattern of the P-IV, peaks are indicated 
by dots and G represents diffraction by a metal gasket. 



side figure is the projection onto the ac plane. It shows 
the distortion of the simulation cell from the cubic into 
the cell with an angle 97.76°. When we look the lat- 
tice from a side, we find a zigzag modulation of the ac 
plane along the 6-axis with displacement in the direction 
of [110] as is shown in the ABAB.. in the right hand side 
figure. This is an important feature. The ABAB- - - mod- 
ulation pattern is crucial for the stability of the distortion 
of the angle /3. When we removed the zigzag modulation 
pattern and performed the simulation for the relaxation 
of the structure, we observed that the structure returned 
to the initial sc structure. 

The structure obtained by the metadynamics is a mod- 
ulation pattern with period consisting of two planes. Our 
simulation, however, was performed using the system 
with 8-atoms in a simulation cell, and with the periodic 
boundary condition. Hence there remains a question that 
the small simulation cell may limit the modulation period 
to a shorter one. In fact, the X-ray powder pattern of the 
zigzag modulated structure and that of the experimental 



one show some discrepancies. To answer this question, 
we tried some more first-principles calculations for the 
refinement of the structure. 

We extended our study to two more structures which 
have commensurate modulations: an ABAC- - - and an 
ABCBADED- - -. The modulation period of the ABAC- - - 
and that of the ABCBADED- - - are twice and four times 
as long as that of the ABAB- - -, respectively. These 
structures are denoted as M4 and M8 and shown in Fig. 
13 The pattern Ml is a non-modulated structure, and 
M2 is the structure obtained by the metadynamics sim- 
ulation. We calculated the total energies and the X-ray 
powder patterns. 

To compare the total energy among these structures, 
we used the same unit cell with the lattice vectors, a/2, 46 
and c/2, where a, b and c are those of the simulation cell 
obtained by the metadynamics simulation. The 8 atoms 
were located only along the b axis in the unit cell with 
displacements corresponding to the modulation pattern 
(Fig. 13) . This choice of the unit cell avoids the numerical 
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errors coming from the use of the different size of the 
unit cells. For the k-space integration, we used 16 x 4 
X 16 mesh points in the first Brillouin zone. Amplitude 
of modulation for each structure was optimized by the 
relaxation of the atomic positions. 

Calculated total energies are listed in Tabled All mod- 
ulated structures, M2, M4, and M8 have lower energy 
per atom than the unmodulated structure, Ml, which 
shows any modulation periods from M2 to M8 are more 
favorable than the Ml. Among the above three modu- 
lated structures, the energy of the M4 is the lowest in 
our study of the commensurate approximation. Though 
the enthalpy of the M4 is very close to that of the sc, it 
is in fact lower than that of the sc at 120GPa, according 
to our results. If we plot the total energy as a function 
of the modulation period and optimize the period as was 
done by Ehlers et al. |l2l| . an incommcnsurately modu- 
lated structure is expected. 

In figure 01 we compare the X-ray powder patterns of 
our structures with that of the experimental one. The ex- 
perimental pattern of the P-IV (top figure) was obtained 
from Akahama et al.'s Q. It shows the feature that the 
splits of the strongest peak at 26 = 13° and the three 
peaks in the range from 29 = 17° to 19° are observed. 
In our M2 structure, these features are missing and un- 
necessary peaks exist at 29 = 15° and 25°. However, 
this disagreement is much improved with increase of the 
modulation period. Intensities of the unnecessary peak at 
29 = 15° in the M2 structure, which was brought about 
by the zigzag modulation, and another unnecessary one 
at 29 ~ 25° in the same M2 structure, which appeared 
owing to the distortion of the simulation cell from the cu- 
bic, are decreased in the M4 and M8 structures and the 
split of the strongest peak appears also in the M4 and 
MB structure. About the intensity of the three peaks in 
the range from 29 = 17° to 19°, the X-ray powder pat- 
tern of the M4 shows most improved agreement with the 
experimental one among the 4 patters studied. From the 
comparison of the total energies and the X-ray powder 
patterns, we conclude that the modulation period is close 
to that of the M4 structure in the P-IV. 

In this study, we explored the structure of phosphorus 
in the phase IV using the first-principles metadynamics 
simulation and identified the new structure. This struc- 
ture is the monoclinic with the modulation pattern. Fur- 
thermore we found the refined structure showing the best 
agreement with the X-ray powder pattern. Although we 
have not fully studied the possibility of the incommen- 
surate modulation, we conclude that phosphorus takes 
the modulated structure with possible incommensurate 
modulation. The simple idea of Akahama et al. stat- 
ing that the structure of the P-IV may be on the path 
from the sc to sh via monoclinic distortion along the [110] 
direction is partly supported because the unmodulated 
structure is of space group P2 and of one atom per unit 



cell. The structure of the P-IV is not so complicated as 
one suggested by and Ehlers et al 0|. The modulation 
stabilizes the monoclinic distortion of the lattice. It is 
highly probable that the Vb, VIb, and Vllb group ele- 
ments commonly show modulated structures in a narrow 
pressure range when they undergo the pressure induced 
structural transition between simple stable structures. 
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